import matplotlib.pyplot as plt
import numpy as np

# Define the logistic growth model as a function
def logistic_growth(P, r, K):
    return r * P * (1 - P / K)


# Parameters for the logistic model
r_value = 1.0  # Growth rate
K_value = 10.0  # Carrying capacity

# Define the range of P values
P_values = np.linspace(0, 12, 500)

# Compute the rate of change dP/dt for each P
dP_dt = logistic_growth(P_values, r_value, K_value)

plt.figure(figsize=(8, 6))
for r_value in np.linspace(0.5, 2.0, 5):
    dP_dt = logistic_growth(P_values, r_value, K_value)
    plt.plot(P_values, dP_dt, label=f"r = {r_value:.1f}")  

#plt.plot(P_values, dP_dt, label=r"$\dot{P} = rP(1 - \frac{P}{K})$", color="blue")
plt.axhline(0, color="black", linestyle="--", linewidth=0.8)  # Zero growth line
plt.axvline(K_value, color="red", linestyle="--", label="Carrying Capacity (K)")
plt.axvline(0, color="red", linestyle="--", label="Extinction (P = 0)")
plt.title("Phase Diagram of Logistic Growth")
plt.xlabel("Population (P)")
plt.ylabel(r"Rate of Change ($\dot{P}$)")
plt.grid(True)
plt.legend()
plt.ylim(-0.1,K_value+0.1)
plt.xlim(-0.1,K_value+0.1)
plt.show()

# Define the function for exponential growth
def exponential_growth(P0, r, t_max, num_points=500):
    t_values = np.linspace(0, t_max, num_points)
    P_values = P0 * np.exp(r * t_values)
    return t_values, P_values

# Define the function for logistic growth
def logistic_growth(P0, r, K, t_max, num_points=500):
    t_values = np.linspace(0, t_max, num_points)
    P_values = K / (1 + (K / P0 - 1) * np.exp(-r * t_values))
    return t_values, P_values

# Define the function for logistic growth with mortality
def logistic_with_mortality(P0, r, K, s, t_max, num_points=500):
    t_values = np.linspace(0, t_max, num_points)
    P_values = []
    P = P0
    dt = t_max / num_points
    for _ in t_values:
        dP_dt = r * P * (1 - P / K) - s * P
        P += dP_dt * dt
        if P < 0:  # Prevent negative population
            P = 0
        P_values.append(P)
    return t_values, np.array(P_values)

# Now recompute and plot all three models
P0 = 10  # Initial population
r = 0.5  # Growth rate
K = 100  # Carrying capacity
s = 0.1  # Mortality rate
t_max = 20  # Maximum time for simulation

# Compute population dynamics
t_exp, P_exp = exponential_growth(P0, r, t_max)
t_log, P_log = logistic_growth(P0, r, K, t_max)
t_log_mort, P_log_mort = logistic_with_mortality(P0, r, K, s, t_max)

# Plot all three models in one panel
plt.figure(figsize=(8, 4),dpi=150)

# Exponential growth
plt.plot(t_exp, P_exp, label="Exponential Growth", color="blue", linestyle="-")

# Logistic growth
plt.plot(t_log, P_log, label="Logistic Growth", color="green", linestyle="--")

# Logistic growth with mortality
plt.plot(t_log_mort, P_log_mort, label="Logistic Growth with Mortality", color="red", linestyle=":")

# Add titles and labels
plt.title("Comparison of Population Growth Models")
plt.xlabel("Time (t)")
plt.ylabel("Population (P)")
plt.ylim(0,2*K)
plt.grid(True)
plt.legend()
plt.xlim(0,20)
# Show the plot
plt.show()

# Import sympy for symbolic computation
import sympy as sp

# Define variables for the Lotka-Volterra system
x, y = sp.symbols('x y')
alpha, beta, gamma, delta = sp.symbols('alpha beta gamma delta')

# Define the equations for the Lotka-Volterra system
dx_dt = alpha * x - beta * x * y  # Prey dynamics
dy_dt = delta * x * y - gamma * y  # Predator dynamics

# Define the system as a matrix
F = sp.Matrix([dx_dt, dy_dt])

# Define the state vector
state_vector = sp.Matrix([x, y])

# Compute the Jacobian matrix
J = F.jacobian(state_vector)

# Display the Jacobian matrix
J

# Compute Jacobian at the first fixed point (0, 0)
J_fixed_1 = J.subs({x: 0, y: 0})

# Compute Jacobian at the second fixed point (gamma/delta, alpha/beta)
J_fixed_2 = J.subs({x: gamma / delta, y: alpha / beta})

# Display the Jacobian matrices for both fixed points
display(J_fixed_1)
display(J_fixed_2)

# Compute the eigenvalues for the Jacobian at the second fixed point
eigenvalues_fixed_1 = J_fixed_1.eigenvals()

# Display the eigenvalues
display(eigenvalues_fixed_1)

# Compute the eigenvalues for the Jacobian at the second fixed point
eigenvalues_fixed_2 = J_fixed_2.eigenvals()

# Display the eigenvalues
display(eigenvalues_fixed_2)

{alpha: 1, -gamma: 1}

{-sqrt(-alpha*gamma): 1, sqrt(-alpha*gamma): 1}

import pandas as pd

# Define the symbolic expressions for the eigenvalues of the Lotka-Volterra system
x, y = sp.symbols('x y')
alpha, beta, gamma, delta = sp.symbols('alpha beta gamma delta')

# Jacobian matrix for Lotka-Volterra at the second fixed point (gamma/delta, alpha/beta)
J_lv_fixed_2 = sp.Matrix([
    [0, -beta * gamma / delta],
    [delta * alpha / beta, 0]
])

# Calculate eigenvalues of the Jacobian at the second fixed point
eigenvalues_lv = J_lv_fixed_2.eigenvals()

# Convert eigenvalues into a numerical framework for plotting
alpha_vals = np.linspace(0, 2, 111)  # Range of alpha values
gamma_val = 1.0  # Fix gamma for this analysis

# Lists to store real and imaginary parts

x_v = np.linspace(0, 2, 100)
df = {"alfa":[], "lambda_1":[], "lambda_2":[]}

l1r = []
l1i = []

l2r = []
l2i = []

for alpha_val in alpha_vals:
    # Substitute alpha and gamma into eigenvalues
    ev_numeric = [ev.subs({alpha: alpha_val, gamma: gamma_val}) for ev in eigenvalues_lv.keys()]
    
    # Real and imaginary parts of each eigenvalue
    lambda_1_real = sp.re(ev_numeric[0])
    lambda_1_imag = sp.im(ev_numeric[0])

    lambda_2_real = sp.re(ev_numeric[1])
    lambda_2_imag = sp.im(ev_numeric[1])

    l1r.append(lambda_1_real)
    l1i.append(lambda_1_imag)

    l2r.append(lambda_2_real)
    l2i.append(lambda_2_imag)

    df["alfa"].append(alpha_val)
    
    sign_of_impart = "+" if lambda_1_imag >= 0 else "-"
    df["lambda_1"].append(f"{lambda_1_real:.2f} {sign_of_impart} {abs(lambda_1_imag):.2f}i")
    sign_of_impart = "+" if lambda_2_imag >= 0 else "-"
    df["lambda_2"].append(f"{lambda_2_real:.2f} {sign_of_impart} {abs(lambda_2_imag):.2f}i")

df = pd.DataFrame(df)
display(df)

# Plot real and imaginary parts of eigenvalues
plt.figure(figsize=(10, 8))

# Real parts
plt.plot(alpha_vals, l1r, label='$Re(\\lambda_1)$', color="navy")
plt.plot(alpha_vals, l1i, label='$Im(\\lambda_1)$', color="cyan")

# Imaginary parts

plt.plot(alpha_vals, l2r, label='$Re(\\lambda_2)$', color="darkred")
plt.plot(alpha_vals, l2i, label='$Im(\\lambda_2)$', color="red")

# Add labels, legend, and grid
plt.axhline(0, color='black', linestyle='dashed', linewidth=0.8)
plt.xlabel('Alpha ($\\alpha$)')
plt.ylabel('Eigenvalue Components')
plt.title('Eigenvalue Components for Lotka-Volterra System')
plt.legend()
plt.grid(True)

plt.show()

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def lotka_volterra_ivp(t, state, alpha, beta, delta, gamma):
    x, y = state
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# Initial conditions
x0 = 10
y0 = 5

# Parameters
alpha = 0.1
beta = 0.4
delta = 0.1
gamma = 0.4
parameters = (alpha, beta, delta, gamma)

# Time span and evaluation points
t_span = (0, 100)
t_eval = np.linspace(t_span[0], t_span[1], 500)

# Solve the IVP
solution = solve_ivp(lotka_volterra_ivp, t_span, [x0, y0], args=parameters, t_eval=t_eval)

# Extract the solution
time = solution.t
x, y = solution.y

# Plot the solution
plt.figure(figsize=(8, 6))
plt.plot(time, x, label='Prey (x)', color="blue")
plt.plot(time, y, label='Predator (y)', color="red")
plt.xlabel('Time')
plt.ylabel('Population')
plt.title('Lotka-Volterra Model')
plt.legend()
plt.grid(True)
plt.show()

# Parameters for the system
alpha_val = 1.0
beta_val = 0.5
gamma_val = 1.0
delta_val = 0.5

# Define the Lotka-Volterra system
def lotka_volterra_ivp(t, state, alpha, beta, delta, gamma):
    x, y = state
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]

# Initial conditions for trajectories
initial_conditions = [
    (1, 1),
    (1, 2),
    (1, 3),
    (1, 4)
]

# Time range for simulation
t_values = np.linspace(0, 50, 1000)

# Create a grid for the vector field
x_values = np.linspace(0, 6, 20)
y_values = np.linspace(0, 6, 20)
X, Y = np.meshgrid(x_values, y_values)

# Calculate the vector field
dX = alpha_val * X - beta_val * X * Y
dY = delta_val * X * Y - gamma_val * Y

# Plot the phase diagram with trajectories
plt.figure(figsize=(10, 8))

# Plot the vector field
plt.quiver(X, Y, dX, dY, angles='xy', scale_units='xy', scale=50, color="gray") 

# Plot trajectories for different initial conditions using solve_ivp
for x0, y0 in initial_conditions:
    # Solve the Lotka-Volterra equations
    sol = solve_ivp(
        lotka_volterra_ivp, 
        [0, 50], 
        [x0, y0], 
        args=(alpha_val, beta_val, delta_val, gamma_val), 
        t_eval=t_values
    )
    plt.plot(sol.y[0], sol.y[1], label=f"Trajectory (x0={x0}, y0={y0})")

# Mark the fixed points
fixed_points = [(0, 0), (gamma_val / delta_val, alpha_val / beta_val)]
for fx, fy in fixed_points:
    plt.plot(fx, fy, 'ro', markersize=8, label=f"Fixed Point ({fx:.2f}, {fy:.2f})")

# Add labels, legend, and grid
plt.title("Phase Diagram of Lotka-Volterra System")
plt.xlabel("Prey Population (x)")
plt.ylabel("Predator Population (y)")
plt.grid(True)
plt.legend()
plt.xlim([0, 6])
plt.ylim([0, 6])

# Show the plot
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import convolve

# Parameters for the Lotka-Volterra model
alpha, beta, gamma, delta = 1.0, 0.1, 1.5, 0.075
Dx, Dy = 0.01, 0.01  # Diffusion coefficients
L = 100  # Grid size
dx = 1.0  # Spatial resolution
dt = 0.1  # Time step
T = 100  # Total simulation time
steps = int(T / dt)

# Initialize prey (x) and predator (y) populations
x = np.random.uniform(0.5, 1.5, (L, L))
y = np.random.uniform(0.5, 1.5, (L, L))

# Laplacian kernel for diffusion
laplacian_kernel = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]]) / dx**2

# Simulation loop
for step in range(steps):
    # Compute reaction terms
    x_react = alpha * x - beta * x * y
    y_react = delta * x * y - gamma * y

    # Compute diffusion terms
    x_diff = Dx * convolve(x, laplacian_kernel, mode='wrap')
    y_diff = Dy * convolve(y, laplacian_kernel, mode='wrap')

    # Update populations
    x += (x_react + x_diff) * dt
    y += (y_react + y_diff) * dt

    # Ensure non-negative populations
    x = np.maximum(x, 0)
    y = np.maximum(y, 0)

    # Visualization at specific intervals
    if step % 500 == 0 or step == steps - 1:
        plt.figure(figsize=(12, 6))
        plt.subplot(1, 2, 1)
        plt.title(f"Prey Population (Step {step})")
        plt.imshow(x, cmap='viridis', extent=[0, L, 0, L])
        plt.colorbar()
        plt.subplot(1, 2, 2)
        plt.title(f"Predator Population (Step {step})")
        plt.imshow(y, cmap='plasma', extent=[0, L, 0, L])
        plt.colorbar()
        plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import convolve
import matplotlib.animation as animation
from IPython.display import HTML

# Parameters
alpha, beta, gamma, delta = 1.0, 0.1, 1.5, 0.075
Dx, Dy = 0.01, 0.01
L = 100  # Grid size
dx = 1.0
dt = 0.1

# Initialize populations
x = np.random.uniform(0.5, 1.5, (L, L))
y = np.random.uniform(0.5, 1.5, (L, L))

# Laplacian kernel
laplacian_kernel = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]]) / dx**2

# Create figure
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
im1 = ax1.imshow(x, cmap='viridis', extent=[0, L, 0, L], animated=True)
im2 = ax2.imshow(y, cmap='plasma', extent=[0, L, 0, L], animated=True)
ax1.set_title("Prey Population")
ax2.set_title("Predator Population")

def update(frame):
    global x, y
    # Compute reaction terms
    x_react = alpha * x - beta * x * y
    y_react = delta * x * y - gamma * y

    # Compute diffusion terms
    x_diff = Dx * convolve(x, laplacian_kernel, mode='wrap')
    y_diff = Dy * convolve(y, laplacian_kernel, mode='wrap')

    # Update populations
    x[:] += (x_react + x_diff) * dt
    y[:] += (y_react + y_diff) * dt

    # Ensure non-negative populations
    x[:] = np.maximum(x, 0)
    y[:] = np.maximum(y, 0)

    # Update plot
    im1.set_array(x)
    im2.set_array(y)
    return im1, im2

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=200, interval=50, blit=True)
ani.save('lotka_volterra_animation.mp4', writer='ffmpeg', fps=30)
# Convert the animation to HTML and display
HTML_animation = HTML(ani.to_html5_video())
HTML_animation

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import convolve
import matplotlib.animation as animation
from IPython.display import HTML
plt.close()
# Parameters
alpha, beta, gamma, delta = 1.0, 0.1, 1.5, 0.075
Dx, Dy = 0.01, 0.01
L = 100  # Grid size
dx = 1.0
dt = 0.1
steps = 500  # Number of time steps for animation

# Initialize populations
x = np.random.uniform(0, 1, (L, L))
y = np.random.uniform(0, 1, (L, L))

# Laplacian kernel
laplacian_kernel = np.array([[0, 1, 0], [1, -4, 1], [0, 1, 0]]) / dx**2

# Time series data
time = []
prey_totals = []
predator_totals = []

# Create figure for animation
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Matrix plot
im1 = ax1.imshow(x - y, cmap='coolwarm', extent=[0, L, 0, L], animated=True)
ax1.set_title("Matrix Plot (Prey-Predator)")
ax1.set_xlabel("Space")
ax1.set_ylabel("Space")
cbar = plt.colorbar(im1, ax=ax1)
cbar.set_label("Population Difference")

# Time series plot
ax2.set_title("Time Series (Prey and Predator)")
ax2.set_xlim(0, steps * dt)
ax2.set_ylim(0, 1.5 * np.max([x.sum(), y.sum()]))
ax2.set_xlabel("Time")
ax2.set_ylabel("Population Sum")
prey_line, = ax2.plot([], [], label="Prey (Sum)", color="green")
pred_line, = ax2.plot([], [], label="Predator (Sum)", color="red")
ax2.legend()

def update(frame):
    global x, y, time, prey_totals, predator_totals

    # Compute reaction terms
    x_react = alpha * x - beta * x * y
    y_react = delta * x * y - gamma * y

    # Compute diffusion terms
    x_diff = Dx * convolve(x, laplacian_kernel, mode='wrap')
    y_diff = Dy * convolve(y, laplacian_kernel, mode='wrap')

    # Update populations
    x[:] += (x_react + x_diff) * dt
    y[:] += (y_react + y_diff) * dt

    # Ensure non-negative populations
    x[:] = np.maximum(x, 0)
    y[:] = np.maximum(y, 0)

    # Update time series data
    time.append(frame * dt)
    prey_totals.append(x.sum())
    predator_totals.append(y.sum())

    # Update matrix plot
    im1.set_array(x - y)

    # Update time series plot
    prey_line.set_data(time, prey_totals)
    pred_line.set_data(time, predator_totals)

    # Set x-axis limits only if there are multiple time points
    if len(time) > 1:
        ax2.set_xlim(0, time[-1])
    ax2.set_ylim(0, 1.5 * max(max(prey_totals), max(predator_totals)))
    return im1, prey_line, pred_line

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=steps, interval=50, blit=True)
ani.save('lotka_volterra_animation_v2.mp4', writer='ffmpeg', fps=30)
# Convert the animation to HTML for Jupyter Notebook
HTML_animation = HTML(ani.to_html5_video())
HTML_animation

C:\Users\Rene\AppData\Local\Temp\ipykernel_7836\622603565.py:81: UserWarning: Attempting to set identical low and high xlims makes transformation singular; automatically expanding.
  ax2.set_xlim(0, time[-1])
C:\Users\Rene\AppData\Local\Temp\ipykernel_7836\622603565.py:81: UserWarning: Attempting to set identical low and high xlims makes transformation singular; automatically expanding.
  ax2.set_xlim(0, time[-1])

	alfa	lambda_1	lambda_2
0	0.000000	0.00 + 0.00i	0.00 + 0.00i
1	0.018182	0.00 - 0.13i	0.00 + 0.13i
2	0.036364	0.00 - 0.19i	0.00 + 0.19i
3	0.054545	0.00 - 0.23i	0.00 + 0.23i
4	0.072727	0.00 - 0.27i	0.00 + 0.27i
...	...	...	...
106	1.927273	0.00 - 1.39i	0.00 + 1.39i
107	1.945455	0.00 - 1.39i	0.00 + 1.39i
108	1.963636	0.00 - 1.40i	0.00 + 1.40i
109	1.981818	0.00 - 1.41i	0.00 + 1.41i
110	2.000000	0.00 - 1.41i	0.00 + 1.41i

Modeliranje rasti populacije ¶

Omejeni viri ¶

Smrtnost ¶

Stabilnost fiksnih točk ¶

Kaj smo vse upoštevali ¶

Lotka-Voltera ¶

Fiksne točke ¶

Stabilna analiza ¶

Solve ode ¶

Kaj smo vse upoštevali ¶

Prostorska interakcija v Lotka-Volterra modelu ¶

Ključne komponente prostorskega modela ¶

Prostorski Lotka-Volterra model ¶

Simulacija prostorske interakcije ¶

Kako začeti simulacijo?¶

Kaj smo vse upoštevali ¶